#clean up
rm(list=ls())

#set seed
set.seed(715130425)

#load libraries
library(akima)
library(geoR)
library(scatterplot3d)  
library(ggmap)
library(lattice)
library(foreign)
library(maps)
library(spdep)
library(doParallel)
library(foreach)

#options
options(scipen=8)

#Register Number of Processors for Parallel Computing
#detectCores() 
registerDoParallel(7)
getDoParWorkers()

#Set to working directory. Adjust the command below to reflect YOUR working directory.
#setwd('/Volumes/MONOGAN/gillAreal/ccesCode/')
#setwd('/Users/jamie/Documents/gillAreal/ccesCode/')
#setwd('/Users/jamie/Desktop/nhgis0003_shape/nhgis0003_shapefile_us_tract_1970/')

#load raw data
combined<-read.dta('cces08reformatted.dta')

#load kriged points
covariates<-read.csv('krigedPoints.csv')

###BAYESIAN MODEL USING BOOTSTRAPPED SAMPLES###
#number of iterations
K<-100
getDoParWorkers()

#start parallel loop
out<-foreach(k=1:K, .combine=list, .inorder=FALSE, .errorhandling='remove') %dopar%{
	
	#create subset of original data
	five.percent<-sample(1:dim(combined)[1],round(.05*dim(combined)[1]))
	subset.data<-combined[five.percent,]

	#draw 5 percent of kriged points
	five.targets<-sample(1:dim(covariates)[1],round(.05*dim(covariates)[1]))
	subset.covariates<-covariates[five.targets,]

	#number of simulated observations
	num.sims<-dim(subset.covariates)[1]

	#estimate model on subset
	ideol.bayes.local <- krige.bayes(coords=cbind(subset.data$eastings,subset.data$northings), data=subset.data$ideology, locations=NULL,
		model = model.control(trend.d =~subset.data$age*subset.data$educ+as.factor(subset.data$race)*subset.data$female+
		as.factor(subset.data$inc14)+subset.data$catholic+subset.data$mormon+subset.data$orthodox+subset.data$jewish+
		subset.data$islam+subset.data$mainline+subset.data$evangelical+subset.data$rural+subset.data$ownership+
		as.factor(subset.data$empstat)+subset.data$eastings*subset.data$northings+I(subset.data$eastings^2)+I(subset.data$northings^2),	cov.model = "gaussian"),
		prior = prior.control(beta.prior="flat", tausq.rel.prior="uniform", tausq.rel.discrete=seq(from=6,to=8,by=.05)))

	#iterative output of posterior
	postSample<-ideol.bayes.local$posterior$sample
	write.csv(postSample,paste('sample',k,'.csv',sep=""),row.names=F)

	#simple kriging by iteration
	#create output matrices
	predictions<-matrix(NA,nrow=1000,ncol=num.sims)
	variances<-matrix(NA,nrow=1000,ncol=num.sims)
	
	#loop through each iteration
	for(i in 1:1000){
		kriges<-krige.conv(coords=cbind(subset.data$eastings,subset.data$northings), data=subset.data$ideology, locations=cbind(subset.covariates$east.K,subset.covariates$north.K),
			krige=krige.control(type.krige='sk', trend.d=~subset.data$age*subset.data$educ+as.factor(subset.data$race)*subset.data$female+
				as.factor(subset.data$inc14)+subset.data$catholic+subset.data$mormon+subset.data$orthodox+subset.data$jewish+
				subset.data$islam+subset.data$mainline+subset.data$evangelical+subset.data$rural+subset.data$ownership+
				as.factor(subset.data$empstat)+subset.data$eastings*subset.data$northings+I(subset.data$eastings^2)+I(subset.data$northings^2),
			trend.l=~subset.covariates$age*subset.covariates$educ+as.factor(subset.covariates$race)*subset.covariates$female+
				as.factor(subset.covariates$ftotinc)+subset.covariates$cath+subset.covariates$mormon+subset.covariates$ortho+subset.covariates$jewish+
				subset.covariates$islam+subset.covariates$main+subset.covariates$evan+subset.covariates$rural+subset.covariates$ownershp+
				as.factor(subset.covariates$empstat)+subset.covariates$east.K*subset.covariates$north.K+I(subset.covariates$east.K^2)+I(subset.covariates$north.K^2),
			beta=as.numeric(postSample[i,1:38]), cov.pars=c(postSample$sigmasq[i],postSample$phi[i]), 
			nugget=I(postSample$sigmasq[i]*postSample$tausq.rel[i]), cov.model='gaussian'))
		predictions[i,]<-kriges$predict
		variances[i,]<-kriges$krige.var
	}
	
	#mean for each location
	predictive.mean<-apply(predictions,2,mean)
	
	#variance for each location
	between<-apply(predictions,2,var)
	within<-apply(variances,2,mean)
	predictive.variance<-within+((num.sims+1)/num.sims)*between
	
	#create matrix of kriged points
	forecasts<-cbind(subset.covariates[,c(1,9:14,24)],predictive.mean,predictive.variance)

	#iterative output of forecasts
	write.csv(forecasts,paste('forecast',k,'.csv',sep=""),row.names=F)

}


